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ABSTRACT 

Numerical weather prediction (NWP) simulations using the Japan Meteorological Agency Nonhydrostatic 
Model (JMA-NHM) are conducted for three precipitation events observed by shipborne or spaceborne 
W-band cloud radars. Spectral bin and single-moment bulk cloud microphysics schemes are employed sep- 
arately for an intercomparative study. A radar product simulator that is compatible with both microphysics 
schemes is developed to enable a direct comparison between simulation and observation with respect to the 
equivalent radar reflectivity factor Ze, Doppler velocity (DV), and path-integrated attenuation (PIA). In 
general, the bin model simulation shows better agreement with the observed data than the bulk model 
simulation. The correction of the terminal fall velocities of snowflakes using those of hail further improves the 
result of the bin model simulation. The results indicate that there are substantial uncertainties in the mass-size 
and size-terminal fall velocity relations of snowflakes or in the calculation of terminal fall velocity of snow 
aloft. For the bulk microphysics, the overestimation of Ze is observed as a result of a significant predominance 
of snow over cloud ice due to substantial deposition growth directly to snow. The DV comparison shows that 
a correction for the fall velocity of hydrometeors considering a change of particle size should be introduced 
even in single-moment bulk cloud microphysics. 


1. Introduction 

Active remote sensing using radar/lidar is useful for 
observing the degree of density of clouds and precip- 
itation. W-band cloud radars for frequencies ranging from 
75 to 111 GHz can observe small particles that are sev- 
eral micrometers in diameter, such as cloud droplets, with 
strong Rayleigh scattering; these particles cannot be 
detected by C-band or X-band precipitation radars. A 
frequency of approximately 94 GHz in W-band radar 
measurements is commonly used to monitor clouds and 
precipitation because of the small attenuation by atmo- 
spheric gas absorption at this frequency (Lhermitte 1987). 
Many observations using ground-based, shipborne, or 


Corresponding author address: Takamichi Iguchi, NASA Goddard 
Space Flight Center, Greenbelt, MD 20771. 

E-mail: takamichi.iguchi@nasa.gov 

DOI: 10.1175/JAS-D-11-0213.1 


airborne W-band radars have been reported by prior 
studies (e.g., Clothiaux et al. 1995). Furthermore, Cloud- 
Sat , equipped with a 94-GHz Cloud Profiling Radar 
(CPR), was launched in April 2006, and its data continue 
to be compiled (Stephens et al. 2002, 2008). Global ob- 
servation by A-Train satellites, including CloudSat , pro- 
vides significant information about the structures of clouds 
and precipitation. Another project, the EarthCARE mis- 
sion, is currently being prepared; both 94-GHz Doppler 
radar and lidar will be loaded on the satellite. It is there- 
fore timely to make a study using datasets from cloud 
radar observations. 

The validation of cloud modeling is a way to make use 
of the measurements of cloud radars. The majority of 
prior studies have used radar-signal simulation packages, 
such as QuickBeam (Haynes et al. 2007), to calculate an 
equivalent radar reflectivity factor Ze from prognostic 
variables, which corresponded to the actually measured 
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Ze; this approach enables a direct signal-based com- 
parison. For example, the global output of a multiscale 
modeling framework (MMF) was compared with Cloud- 
Sat CPR observations (Marchand et al. 2009). In a simu- 
lation of the Madden-Julian oscillation (MJO) using a 
global cloud-resolving model (GCRM) with a horizontal 
grid size of several kilometers, measurements of tropical 
clouds and precipitation were compared with datasets 
from the CloudSat CPR and Tropical Rainfall Measuring 
Mission Precipitation Radar (TRMM PR) (Miura et al. 
2007; Masunaga et al. 2008) as well as remote sensing data 
from the Cloud-Aerosol Lidar and Infrared Pathfinder 
Satellite Observation (CALIPS 0)1 CloudSat satellites 
(Inoue et al. 2010; Satoh et al. 2010). The global dis- 
tributions and structures of clouds and precipitation 
in the Met Office global model were evaluated using 
CloudSat CPR data (Bodas-Salcedo et al. 2008). The 
vertical structures of midlatitude and tropical clouds 
observed using shipborne radar and lidar were com- 
pared with those simulated by the Center for Climate 
System Research-National Institute for Environmental 
Studies (CCSR-NIES) atmospheric general circulation 
model (AGCM) (Okamoto et al. 2007, 2008). The radar 
reflectivity and Doppler velocity of arctic mixed-phase 
clouds simulated by a three-dimensional cloud-resolving 
model coupled with spectral bin microphysics were com- 
pared with measurements by a millimeter cloud radar 
(Fan et al. 2009). 

In this study, we conducted numerical weather predic- 
tion (NWP) simulations using the Japan Meteorological 
Agency Nonhydrostatic Model (JMA-NHM) and com- 
pared the results with observations by shipborne or 
spaceborne cloud radars for three specific midlatitude 
cases. Spectral-bin-microphysics and single-moment-bulk 
schemes for clouds were separately used for an inter- 
comparison study. The corresponding variables in the 
measurement products, such as Ze, Doppler velocity 
(DV), and path-integrated attenuation (PI A), were cal- 
culated using a newly developed radar product simulator 
that had been optimized for the bin scheme. This ap- 
proach may allow a signal-based comparison between 
simulated data and observed data (Masunaga et al. 2010). 
The assumptions in the simulator are highly consistent 
with those made in the model microphysical schemes. 
Our simulator directly utilizes particle size distributions 
(PSDs) of hydrometeors resulting from the bin model 
simulations. Built-in PSDs, assumed originally in the bulk 
microphysics, are introduced into the simulator when 
applied to the bulk model simulations. 

In the first case, data of a 95-GHz Doppler radar on 
a research vessel during a cruise over the Pacific Ocean 
near Japan in May 2001 were used. The second and third 
cases involved cloud systems that occurred near Japan in 


November 2006 that were recorded by CloudSat global 
observation using the 94-GHz CPR. Comparative anal- 
yses of the three cases may provide a better un- 
derstanding of the differences between the simulated and 
observed data, despite the differences in the weather 
conditions and locations, because the zenith and nadir 
remote sensing data are complementary to each other 
with differing attenuation profiles. 

The methodology is described in section 2 and includes 
a description of the NWP simulation and an explanation 
of the radar product simulator and radar observations. 
The comparative analysis between the numerical exper- 
iments and the observed data is described in section 3. 
The summary and conclusions are presented in section 4. 

2. Methodology 

a. Cloud radar observations 

1) MR01/K02 CRUISE OF THE RESEARCH VESSEL 
Mirai 

The MR01/K02 research cruise was conducted with 
a shipborne 95-GHz Doppler radar and lidar with a 
zenith looking over the northwest Pacific near Japan 
from 14 to 28 May 2001 (Sugimoto et al. 2002; Okamoto 
et al. 2007). The observed data provided Ze, linear de- 
polarization ratio (LDR), DV, and the correlation co- 
efficient between the horizontal and vertical polarization 
signals; the technical description can be found in Horie 
et al. (2000). The calibration accuracy of the radar signal is 
approximately ±1 dB (Okamoto 2002; H. Kuroiwa 2001, 
personal communication). The vertical range of the data 
reaches a height of approximately 12 km with intervals of 
82.5 m and a time interval of 1 min. This study targets 22 
and 23 May 2001, when the vessel stayed at approximately 
34°N, 146°E, and a low pressure system typical of the 
midlatitudes with associated fronts approached from the 
west and passed over the vessel (Fig. la). The radar ob- 
served the vertical structures of the clouds and the pre- 
cipitation attributed to the low pressure and frontal 
system. The passing clouds were altostratus, nimbostra- 
tus, and shallow cumulus. 

2) CloudSat CPR observation 

CloudSat is a polar-orbiting satellite with a mean 
equatorial altitude of 705 km in the formation of the five 
satellites known as the A-Train constellation (Stephens 
et al. 2002, 2008). The vertical interval of the CPR data is 
240 m, and the horizontal cross-track and along-track 
resolutions are 1.4 and 1.8 km, respectively. The Ze in 
the CloudSat Geometric Profile (2B-GEOPROF) prod- 
uct (Mace et al. 2007; Marchand et al. 2008; Stephens 
et al. 2008; Tanelli et al. 2008) and the PIA due to 
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Fig. 1. Surface weather charts at (a) 0000 UTC 23 May 2001, (b) 1200 UTC 11 Nov 2001, and (c) 1200 UTC 14 Nov 2001 with overlaid 
black squares corresponding to the simulation domains. The dots and lines denote the coordinates of Mirai and the footprints of CloudSat, 
respectively. (The surface weather charts are provided by the Japan Meteorological Agency.) 


hydrometeors in the Precipitation Column (2C-PRECIP- 
COLUMN) product (Haynes et al. 2009) were provided 
from the CloudSat Data Processing Center (DPC). 

We analyzed two events, when CloudSat passed over 
clouds near Japan on 11 and 14 November 2006 under 
dissimilar weather conditions and in different regions. 
The corresponding CloudSat track granules were num- 
bered 2873 and 2917 by the CloudSat DPC. In the first 
case, the satellite passed over the northeastern part of 
Japan at approximately 1638 UTC 11 November 2006 
from 34.6°N, 141.7°E to 46.1°N, 145.7°E (the latitude and 
longitude boundaries follow those of the QuickLook 
segment provided by the DPC). The CPR observed very 
thick nimbostratus and precipitation located near the 
center of a low pressure system in the mature stage (Fig. 
lb). A uniform vertical structure of reflectivity extended 
from 39° to 46°N along the footprint of the satellite. In the 
second case, the satellite passed over the Japan Sea at 
roughly 1708 UTC 14 November 2006 from 34.6°N, 
134.0°E to 46.1°N, 137.6°E. The CPR observed the 
structures of shallow convective clouds extending behind 
a cold front over the sea surface (Fig. lc). 

b. Description of JMA-NHM 

The operational version of JMA-NHM (Saito et al. 
2006) is the main model framework for three-dimensional 
NWP simulations. The basic governing equations are fully 
compressive nonhydrostatic equations. The time-splitting 
horizontally explicit and vertically implicit (HE-VI) 
scheme is used to inhibit the inflation of sound waves. The 
Arakawa-C and Lorenz grid structures are employed 
in the horizontal and vertical grid discretizations, re- 
spectively. For points that differed from the original JMA- 
NHM, the broadband radiative transfer code “mstrn-x” 
(Sekiguchi and Nakajima 2008) was employed to calculate 


atmospheric radiation. The subgrid convective parame- 
terization scheme in the JMA-NHM was not employed to 
prevent a conflict with the spectral bin scheme, which is 
described below. 

The spectral bin microphysics scheme, which is based 
on the module package of the Hebrew University Cloud 
Model (HUCM; e.g., Khain et al. 2000), is integrated 
into the JMA-NHM (Iguchi et al. 2008). Hydrometeors 
are categorized into one water class and six ice classes 
[i.e., water droplets, ice crystals (plate, column, dendrite), 
snowflakes, graupel, and hail]. Ice crystals are defined as 
primary ice particles that have not coagulated with any 
other hydrometeor particles; snowflakes are aggregates of 
ice crystals. Graupel and hail are rimed particles that 
primarily arise from the combination of supercooled 
droplets and ice hydrometeor particles or from the 
freezing of supercooled droplets. Snowflakes, graupel, 
and hail are assumed to be spheres when calculating their 
microphysics (Khain and Sednev 1995). The discrete 
PSDs of hydrometeors are represented on a grid con- 
taining 33 doubling mass bins 1 covering particles mass 
sizes in a range in which 3.35 X 10 -11 < m < 1.44 X 
10 _1 g (2 < r < 3251 pirn in terms of the radii of droplets 
or melted ice). The scheme calculates nucleation for 
droplets and ice crystals, condensation and deposition 
growths, evaporation, sublimation, droplet freezing, 
melting, and coalescence growth. The ice nucleation rate 
was updated using the equation in Cotton et al. (1986). 

Additional experiments using the bin scheme were 
conducted together: the fall velocities of snowflakes 
were replaced with those of hail. This adjustment was 


1 The sequence of particle mass is a geometric progression with 
common ratio 2. 
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adopted from the latest update of the bin microphysics 
(Khain et al. 2011), in which the fall velocities of snow- 
flakes are variable because of riming effects. This type of 
adjustment for the fall velocity of snow was tested to 
ameliorate the discrepancy between the model simula- 
tions and radar observations in several prior studies 
(Fan et al. 2009; Li et al. 2010). 

The bulk microphysical scheme originally included in 
the JMA-NHM was also integrated to allow a comparative 
analysis. This scheme is a single-moment bulk accounting 
for explicit classes, two of water and three of ice: cloud 
water, rain, cloud ice, snow, and graupel (Lin et al. 1983; 
Ikawa and Saito 1991; Eito and Aonashi 2009). The 
amount of cloud water is calculated using an instantaneous 
saturation adjustment. The autoconversion rate of rain 
from cloud water is calculated by Kessler’s autoconversion 
formula through the collision and coalescence of droplets. 
Cloud ice is defined as pristine ice crystals, and snow is 
defined as snow crystals or aggregates. Graupel is a type of 
rimed ice particle that includes hail. Exponential distri- 
bution functions are assumed to capture the forms of the 
PSDs of rain, snow, and graupel. Monodispersive func- 
tions are hypothesized to characterize cloud water and 
cloud ice. The gravitational sedimentations of particles 
categorized as cloud water and cloud ice are not calculated 
because of their tiny terminal velocities. 2 

c. Setup of the NWP simulations using JMA-NHM 

Three sets of NWP simulations were prepared sepa- 
rately for 22-23 May 2001, 11 November 2006, and 14 
November 2006 (hereinafter, the three cases are referred 
to as M2223, C2873, and C2917). The simulation domains 
for M2223, C2873, and C2917 were centered at 34°N, 
146°E; 42°N, 143°E; and 38°N, 135°E, respectively. The 
simulation domains for all the cases had a horizontal scale 
of 600 km with a grid interval size of 3 km. The vertical 
grid component up to a height of 22 600 m was divided into 
40 layers with intervals increasing from 40 to 1120 m with 
altitude. A time step of 20 s was set; a variable time step 
shorter than this interval was used for the bin microphysics. 

In the case of M2223, nine individual NWP simulations 
with 6 h of integration were conducted for 54 h in total 
from 1800 UTC 21 May to 0000 UTC 24 May 2001. Each 
simulation transferred the specific prognostic variables 
[i.e., potential temperature, mixing ratio of vapor, and 
PSDs of hydrometeor particles and condensation nuclei 
(CN)] to the next 6-h simulation. We continuously sampled 
the prognostic variables in the column over the horizontal 


2 The sedimentation of cloud ice is included in the latest version 
of the JMA-NHM to prevent excessive accumulation in the upper 
layers (Saito et al. 2007). 


grid nearest to the coordinates of the vessel per 1 min 
during 22 and 23 May. In the cases of C2873 and C2917, the 
pairs of two individual simulations with 6 h of integration 
were conducted for 12 h in total from 0600 to 1800 UTC 
on 11 and 14 November 2006, respectively. We sampled 
the prognostic variables in the columns over the near- 
est horizontal grids to the footprints of CloudSat at 
1638 UTC 11 November and at 1708 UTC 14 November. 

To set the initial and lateral boundary conditions of 
the specific prognostic variables (i.e., the two compo- 
nents of horizontal velocities, potential temperature, 
and mixing ratio of the water vapor), the JMA meso- 
analysis dataset (JMA-MANAL) distributed by the 
Japan Meteorological Business Support Center (JMBSC) 
was employed in NWP simulations. This dataset had 
a horizontal grid interval of 10 km, 20 pressure-plane 
levels, and 4 (for 2001) and 8 (for 2006) samples per day. 
The National Centers for Environmental Prediction 
(NCEP) reanalysis data, which had a horizontal resolu- 
tion of T62 Gaussian grid and four samples per day, were 
used to set the sea surface temperature (SST). A one-way 
nesting technique for CN was applied for the bin model 
simulations (Iguchi et al. 2008). The CN fields were cal- 
culated from the spectral radiation transport model for 
aerosol species [the Spectral Radiation Transport Model 
for Aerosol Species (SPRINTARS); e.g., Takemura et al. 
2005] coupled with the CCSR-NIES AGCM (Hasumi 
and Emori 2004), with a horizontal resolution of T106 
Gaussian grid, 20 vertical layers, and 1 sample per day. 

d. Radar product (Ze, Doppler velocity, and PI A) 
simulator 

A radar product simulator was used to convert the 
prognostic variables obtained from the NWP simulations 
to a corresponding product of radar measurement. The 
prognostic variables in the sampled columns are inter- 
polated on the vertical layers of the radar measurements 
before the calculation of the radar simulator. Our sim- 
ulator is based on that developed for comparing ship- 
borne radar/lidar measurements with AGCM simulations 
(Okamoto et al. 2007, 2008). The simulator was modified 
so that it could be plugged into both the bin and bulk 
microphysics of JMA-NHM. The assumptions in the sim- 
ulator (i.e., PSD and particle fall velocities) are consistent 
with those made in the model microphysical schemes. 

The equivalent radar reflectivity factor is calculated 
using the following equation (Okamoto et al. 2003): 


Ze 


A 4 

7t 5 |K| 2 


pr 

max 

Jr . 

min 


dn(r) 

dr 


c bk( r ) dr > 
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where A is the wavelength, r is the radius of particles, n(r) is 
the number concentration of particles with a radius less 
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than r, and C bk is the backscattering cross section given 
from the lookup table that was precalculated on the basis 
of the discrete dipole approximation (DDA) approach. 
Also, K is the dielectric factor estimated from n by \K\ = 

| (n 2 — 1 )/(n 2 + 2)|, where n is the complex refractive in- 
dex. The value of |i£| is set to be 0.828 for the Mirai 
case (Horie et al. 2000; Okamoto et al. 2007) and 0.75 for 
the CloudSat case (Stephens et al. 2008). The different 
K values are derived from their respective equations to 
convert the received power to Ze. Here \K\ is assumed to 
be the same for both water and ice; this assumption is 
based on the definition of Ze for each observation. We use 
the logarithmic form, dBZe = 10 logi 0 Ze, expressed in 
terms of decibels (dB), in the following analysis. 

The backscattering cross sections of all ice hydro- 
meteor particles were calculated using the DDA ap- 
proach regardless of their effective or bulk densities in 
the cloud microphysical schemes. The mass equivalent 
volume with a density of 0.9 g cm -3 was defined for 
each mass bin, and the backscattering cross section was 
computed using DDA for the sphere particle with the 
equivalent volume (Okamoto 2002; Okamoto et al. 
2003; Sato and Okamoto 2006). The sensitivity of the 
result of the radar simulator to spherical and non- 
spherical particle models in the DDA computation is 
briefly discussed in subsection 3a. 

The prognostic PSDs of the hydrometeors in the bin 
model simulation are directly substituted for Eq. (1). In 
the case of the bulk model simulation, the corresponding 
PSDs are calculated using the prognostic mixing ratios 
and built-in PSDs assumed originally in the bulk 
microphysics; then, the PSDs are directly substituted for 
Eq. (1), as in the case of the bin model simulation. The 
total mass concentration is conserved when the PSDs 
are expanded on the bins. The built-in PSDs of rain, 
snow, and graupel are assumed to be exponential dis- 
tributions in the bulk microphysics in the following way: 

N(D) = N 0 exp(-AD), (2) 

where D is the particle diameter, N 0 is the intercept 
parameter (i.e., N r0 = 8.0 X 10 6 m -4 for rain, N s0 = 1.8 X 
10 6 m -4 for snow, and Ngo = 1.1 X 10 6 m -4 for graupel; 
Ikawa and Saito 1991), and A is the slope parameter de- 
termined by the mass mixing ratio with the equations of 
Lin et al. (1983). The particle densities, p r = 1.0 X 
10 3 kg m -3 for rain, p s = 8.4 X 10 1 kg m -3 for snow, and 
p g = 3.0 X 10 2 kg m -3 for graupel (Eito and Aonashi 
2009), are assumed to calculate the slope parameters. 
The PSDs of cloud water and cloud ice are assumed to 
be monodispersive distributions in the bulk micro- 
physics, and their diameters are defined in the follow- 
ing form: 


a) Bin model particle falling velocities 
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Fig. 2. Relationships between the terminal fall velocities at 
1000 hPa and the bulk radii of cloud particles in (a) the bin and (b) 
the bulk microphysics. 
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where q is the mixing ratio and p a is the air density. The 
prescribed number concentrations are N 0c = 1.0 X 
10 8 m -3 for cloud water (Ikawa and Saito 1991) and N 0i = 
2.0 X 10 6 X exp(-0.122T) m -3 for cloud ice (T is the 
temperature in degrees Celsius) (Wilson and Ballard 1999; 
Bodas-Salcedo et al. 2008) (this alternative parameter for 
cloud ice was used because no assumption was made in the 
original). The particle densities are p c = 1.0 X 10 3 kg m -3 
for cloud water and p t = 5.0 X 10 2 kg m -3 for cloud ice. 

The attenuation due to atmospheric hydrometeors for 
Ze is calculated in the following way: 


Ze H = Ze X exp(-2r), (4) 


where r is the optical thickness, that is, the integral of the 
extinction coefficient at the radar wavelength from the 
radar to the target, expressed as (e.g., Kikuchi et al. 2006): 


T 



"( r ) C ext( r ) d r dz , 


(5) 
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a) Observation by ship-borne radar 
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c) Bin model simualtion (rimed snow)+ radar simulator 
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d) Bulk model simulation + radar simulator 
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Fig. 3. Time-height cross sections (THCSs) of the equivalent radar reflectivity factor (dBZe) 
(a) measured by the 95-GHz Doppler radar onboard Mirai, and calculated by the radar product 
simulator applied to the outputs of (b) the bin (control), (c) the bin with the terminal fall 
velocities of snow equalized to those of hail in all size bins (virtually rimed snow), and (d) the 
bulk model simulations from 1200 UTC 22 May to 1200 UTC 23 May 2001. 


where the C ex t values are the extinction cross sections of 
the hydrometeor particles. Also, z r and z t are the heights 
of the radar and the target, respectively; the heights of 
the radars are the surface in a shipborne observation 
case and 705 km aloft in a spaceborne case. The at- 
tenuation resulting from the atmospheric hydrometeor 
components was modeled in the simulator; that is, 
we used radar-measured Ze values without the 


attenuation correction in our comparative analyses. 
The PIA due to hydrometeors is calculated as the two- 
way integrated extinction (Haynes et al. 2009) by 
multiplying r by 2 and 10/ln(10). 

Gas absorption for the typical atmosphere at mid- 
latitude regions is simply parameterized using the 
calculation result in Hogan and Illingworth (1999). The 
total two-way gas attenuation between space and surface 
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a) Observation b) Bin model 
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c) Bin model (rimed snow) 
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Fig. 4. Normalized dBZe-height histograms constructed (a) from the radar measurement and from the simulations 
by (b) the bin (control), (c) the bin (rimed snow), and (d) the bulk model during the same period as Fig. 3. 


is assumed to be 3 dB, and the value from the radar to the 
target A g is formulated in the following form: 

A g = 3 X (1 - e~ zlH )(dB) (shipborne radar), (6) 

A g = 3 X e~ z/H (dB ) (spaceborne radar), (7) 

where H = 2.735 km is the scale height estimated from 
the result in Hogan and Illingworth (1999). 


Only in the case of shipborne radar on the Mirai , the 
radar signal is attenuated by water on the transparent 
cover (radome) of the radar container, which results from 
stagnant precipitation. This attenuation effect was cor- 
rected for the measurement product using the rainfall 
measurement data on the ship; 9 dB was added uniformly 
in the all vertical layers when the measured precipitation 
rate was more than 0.01 mm min -1 following Okamoto 
et al. (2007) (hereinafter, called radome correction). 
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Doppler velocity v d , defined as the sum of reflectivity- 
weighted particle fall velocity and vertical wind velocity, 
is given in the modified form of Matrosov et al. (1994): 

Fmax 

dn(r)/dr • C bk (r)uy(r) dr 

V d = w + )r - , (8) 

u I max 

dn(r)/dr • C bk (r) dr 

^rnin 

where w is the vertical wind velocity and is the terminal 
fall velocities of the hydrometeor particles. In the bin 
microphysics case, Vf is given in the following form: 

v f (r) = -V t (r)^pjp, (9) 

where p is the pressure, p 0 stands for the standard reference 
pressure set at 1000 hPa, and V t is the terminal fall veloc- 
ities of hydrometeor particles at p 0 according to their sizes 
and types (Khain and Sednev 1995) (Fig. 2a). In the case of 
the bulk microphysics, vj- is given by the following equation: 

v f (r) = -a(0.5rf(p aQ lp a y, (10) 

where a, /3, and y are the constants, which can be found 
in Eito and Aonashi (2009). The relationship between 
velocity and size is summarized in Fig. 2b. Similar to the 
calculation of Ze, Eq. (10) with its built-in PSDs is directly 
substituted for Eq. (8) after it has been expanded to a 
spectrum on the size bins. 

3. Results 

a. Mirai shipborne Doppler radar (M2223) 

Figure 3 illustrates the time-height cross sections 
(THCSs) of dBZe, as measured by the 95 GHz radar 
and calculated through the radar simulator with both the 
bin and bulk JMA-NHM NWP simulations for the 24-h 
period from 1200 UTC 22 May to 1200 UTC 23 May 
2001. The first half of 22 May is not included in the 
analysis because few clouds were present during that 
period, and the second half of 23 May is also not in- 
cluded because of the heavy precipitation. The ice and 
water hydrometeor layers were partitioned at a height of 
approximately 4 km. The freezing level was estimated 
from the temperature profile in the mesoanalysis dataset 
and the LDR profile measured by the shipborne radar. 
The upper ice clouds, the lower shallow water clouds, and 
precipitation from the ice clouds were continuously 
observed. This overall comparison shows some spatial or 
temporal mismatches between the measurement and 
simulation because the data are sampled only on the 
particular spot. The total prediction error including the 




Radar sig [dBZe] 

Fig. 5. Vertical distribution of time-averaged dBZe with sampling 
of the range from -40 to 20 dB. Solid, dash-dotted, dash-double 
dotted, and dotted lines denote distributions by the radar measure- 
ment, the bin (control and rimed snow) simulations, and the bulk 
model simulation, respectively, (b) Black (sphe; control), red (br3d), 
blue (cb50), green (cl2d), orange (cl3d), purple (pl2d), and sky blue 
(pl3d) correspond to distributions using the following models: 
spherical, bullet-rosette oriented in three-dimensional space (3D 
bullet-rosette), 50/50 mix of 2D column and 3D bullet-rosette model, 
hexagonal column oriented in horizontal plane (2D column), hex- 
agonal column oriented randomly in three-dimensional plane (3D 
column), and hexagonal plate with 2D orientation (2D plate) and 
random 3D orientation (3D plate) in the DDA computation for 
backscattering cross sections. 

spatial and temporal mismatches should be evaluated in 
the comparison. However, the mismatches are often under 
the control of the dataset to make initial and boundary 
conditions in the NWP simulations rather than the 
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Fig. 6. Vertical distribution of time-averaged mass ratios of each hydrometeor category to the total water contents 
during the same period as Fig. 3. (a) Droplets (solid line), the total of three types of ice crystals (dash-dotted), snow 
(dotted), and the total of graupel and hail (dash-double dotted) in the bin model simulation; (b) the total of cloud 
water and rain (solid line), cloud ice (dash-dotted), snow (dotted), and graupel (dash-double dotted) in the bulk 
model simulation. 


capability of the model physics. To evaluate the perfor- 
mance of the model cloud microphysics, relaxation of the 
error is attempted by obtaining a temporal average of the 
result. A quantitative comparison is shown in the forms of 
the normalized dBZe-height histogram (Fig. 4) and the 
line chart describing the vertical profiles of the 24-h-av- 
eraged dBZe that sampled the range from -40 to 20 dB 
(Fig. 5a). The radome correction for the dBZe measure- 
ment contributes +0.8 dB to the 24-h-averaged dBZe. 

Under the freezing level of approximately 4 km, the 
Ze profiles of the simulations share some common 
characteristics with that of the observed data. The high- 
frequency dBZe bands in the observed data and the bin 
and bulk simulations converge from approximately —10 
to +10 dB near the surface in the histograms, and the 
averaged dBZe shifts to a lower value with altitude. The 
under estimations of dBZe in the bin model simulations 
are seen from 2 to 4 km, with an approximately 5-dB 
difference on average (Fig. 5a). This difference is con- 
sidered to be due to a partial insufficiency of pre- 
cipitation from the ice clouds, which is originally caused 
by a prediction error in the height of ice clouds from 
1200 to 2000 UTC 22 May (Fig. 3). In contrast, over- 
estimations of dBZe are highlighted in the ice cloud 
layers over 4 km in the bin and bulk model simulations. 
The maximum differences in the 24-h-averaged dBZe 
are approximately 5 dB near the height of 5 km for the 
bin model simulation and more than 12 dB for the bulk 


model simulation. This overestimation of dBZe is can- 
celed in the bin model simulation considering the fall 
velocity correction of snow (Figs. 3c, 4c, and 5a). An 
increase in fall velocities led to a decrease in dBZe 
through the decrease in both ice water content (IWC) 
and the mean particle size, because large-size ice hy- 
drometeors were removed more quickly from the layer. 

Figure 5b shows the sensitivity of the 24-h-averaged 
dBZe to different ice particle models in the calculation of 
the backscattering cross section. The following six non- 
spherical particle models were also chosen as 3D particle 
models in the DDA computation: a bullet-rosette 
oriented in three-dimensional space (3D bullet-rosette), 
a hexagonal column oriented in a horizontal plane (2D 
column), a hexagonal column oriented randomly in a 
three-dimensional plane (3D column), a hexagonal plate 
with a 2D orientation (2D plate) and random 3D orien- 
tation (3D plate), and a 50/50 mix of the 2D column and 
the 3D bullet-rosette model (CB50). The details of the 
particle models are given by Sato and Okamoto (2006). 
The backscattering cross sections of the fractal-shaped 
snowflakes at 95 GHz are estimated to be best fitted with 
those of the 2D column model among the spherical and 
the six nonspherical models (Ishimoto 2008; Okamoto 
et al. 2010). Figure 5b shows that the difference in dBZe 
is roughly within the maximum range from —1 to +7 dB, 
and the difference between the spherical model and the 
2D column (considered compatible with fractal snow) 
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a) Observation by ship-borne radar 
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Fig. 7. THCSs of the DVs (a) measured by the 95-GHz Mirai Doppler radar, and calculated 
from the radar product simulator applied to the outputs of (b) the bin (control), (c) the bin (rimed 
snow), and (d) the bulk model simulations from 1200 UTC 22 May to 1200 UTC 23 May 2001. 


model cases is at most 4 dB; the differences depend on the 
particle size and IWC so that the value shifts smaller with 
altitude. This result demonstrates that the utilization of 
the spherical particle model is not a source of the dBZe 
overestimation in the ice cloud layer. 

Figure 6 shows the vertical distribution of the dominant 
hydrometeor types in both the bin and bulk simulations. 
These profiles indicate the comparative characteristics of 
both microphysics and the hydrometer types that caused 


the discrepancies between the observed data and the 
simulations. The overdominance of snow up to the cloud 
top in the bulk model simulation is consistent with other 
results using a similar microphysical framework (Hashimoto 
et al. 2007; Eito and Aonashi 2009). The contoured area 
of the high frequency in the dBZe-height histogram of 
the bulk model simulation over 4 km is small in spread 
(Fig. 4d), and a weak relationship between height and 
dBZe is seen compared with those of both the observed 
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Fig. 8. Normalized DV-height histograms constructed from (a) the radar measurement, (b) the bin (control), (c) the 
bin (rimed snow), and (d) the bulk model simulation. 


data and the bin model simulation. Bodas-Salcedo et al. 
(2008) indicated that this lack of spread is due to an ap- 
plication of the bulk scheme with an assumption of built- 
in PSD to the simulator. In contrast, the corresponding 
spread is sufficiently wide in the histogram of the bin 
model simulation, suggesting that the direct utilization of 
prognostic PSD and the distribution of dominant hydro- 
meteor categories may result in a better dBZe simulation. 
The weak dBZe-height relationship in the bulk model 


simulation is caused by using the constant intercept pa- 
rameter for snow without modeling the dependence on 
temperature, which can be seen from the comparison 
with the results in Bodas-Salcedo et al. (2008). In addi- 
tion, this monotonic Ze structure in the bulk model 
simulation may likely be attributable to the remarkable 
predominance of snow over cloud ice (Fig. 6b). 

An analysis of DV offers an additional perspective to 
evaluate the differences between the observed data and 
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Doppler velocity [m/s] 

Fig. 9. Vertical distribution of the time-averaged Doppler ve- 
locity by radar measurement (solid), bin simulation (control) 
(dash-dotted), bin simulation (rimed snow) (dash-double dotted), 
and bulk model simulation (dotted). 


the simulations (Figs. 7, 8, and 9). Note that the negative 
velocities are downward, and the radar-observed and 
simulated DVs are shown only on vertical areas of the 
dBZe larger than -20 dB. The characteristics of vertical 
structures are separated between the water and ice cloud 
layers, similarly to those of Ze. The DV-height histograms 
of the observation and the bin model simulations share 
some common characteristics: a unimodal distribution is 
dominant over the freezing level of roughly 4 km aloft, 
while a bimodal distribution is seen under the level. The 
small DV mode at approximately -3.5 m s -1 in the bi- 
modal distribution corresponds to the precipitation from 
ice clouds over the freezing level, whereas the large DV 
mode at approximately -1ms -1 corresponds to warm 
clouds whose cloud tops do not reach the freezing level. 
The very weak dominance of the small DV mode in the 
bin model simulation most likely indicates a partial in- 
sufficiency of precipitation from the ice clouds, as in the 
dBZe comparison. 

The bulk model simulation clearly underestimated the 
DV. The velocity decreases with altitude in the water 
and ice cloud layers individually because the low air 
density in the upper atmosphere results in an increase in 
the particle fall velocity [Eq. (10)]. The small particle 
classes such as cloud water and cloud ice, whose fall 
velocity is nearly zero, has little contribution to forming 
DV, probably because these classes are less dominant 
and the monodispersive PSDs with small particle sizes 
cause very small reflectivity. The vertically constant DV 


— ♦ — Observation 

- ■«- • Bin model 

- - -a- - ■ Bin model (only particle) 
— * — Bin model (only air) 

- -* - Bin model - rimed snow 
— • — Bulk model 



Standard deviation of Doppler velocity [m/s] 

Fig. 10. Vertical distribution of the standard deviations of 
Doppler velocity from the time average by radar measurement 
(diamonds), in the bin simulation (control) (squares), only from the 
terminal fall velocities of hydrometeor particles in the bin simu- 
lation (triangles), from vertical wind velocity in the bin model 
simulation (control) (crosses), bin simulation (rimed snow) (as- 
terisks), and bulk model simulation (circles). 


in the bin model simulations suggests that the decrease 
of DV with altitude should be offset by a decrease in the 
fall velocity through a decrease in particle size in the 
cloud microphysics. However, this effect does not work 
properly in the bulk model simulation. 

Even the bin model simulation cannot reproduce the 
observed wide spread over 4 km in the DV-height his- 
togram. A poor reproducibility of the vertical wind ve- 
locity to DV in the simulation could be a source of the 
discrepancy, although the velocity averaged over a long 
period at a given altitude would be expected to converge 
to almost zero in uniform stratiform cloud cases (Orr 
and Kropffi 1999). Figure 10 shows vertical profiles of 
the standard deviation of D V from their averaged values 
in the observation and the simulations. A characteristic 
common to all profiles is a jump at approximately 5 km 
(i.e., the boundary between the ice and water cloud 
layers). This jump is considered to be due to a difference 
in the dependence of the terminal fall velocity on par- 
ticle size between liquid-phase and ice-phase particles 

(Fig. 2). 
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Fig. 11. Normalized dBZe-DV histograms under heights of 4 km constructed from (a) the radar measurement, 
(b) the bin (control), (c) the bin (rimed snow), and (d) the bulk model simulation. 


There is a somewhat constant difference between the 
profiles of the observation and simulations regardless of 
altitude, except over 10 km. This difference could be 
related to the difference of the spatial dispersion of the 
vertical wind velocity between the real atmosphere and 
the model simulations. Although an accurate value of 
the standard deviation contributed by the real wind 
velocity cannot be calculated in this analysis, a value on 
the order of 0.1 ms -1 is feasible compared with the 


retrieval result of the measurement using the same radar 
for tropical cirrus clouds (Sato et al. 2009). In contrast, 
both the standard deviation and the average (not fig- 
ured) of the vertical wind velocity in the model simula- 
tions are approximately less than 0.1 ms -1 at any 
altitude. The difference is probably due to the coarse 
horizontal and vertical resolution of the present NWP 
simulation. We speculate that a horizontal resolution on 
the order of 10 or 1 m is required to reproduce a realistic 
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Fig. 12. As in Fig. 11, but over heights of 5 km. 


dispersion of the vertical wind velocity on a grid-scale 
simulation. Khairoutdinov and Randall (2003) showed 
the sensitivity of vertical velocity variance to the hori- 
zontal resolution in the 28-day simulation during sum- 
mer over Southern Great Plains area. The magnitude of 
variance in the 250-m resolution run is approximately 
twice as much as that in the 2-km run. Although the 
conditions of their simulation are largely different from 
ours, much higher horizontal resolution than 250 m is 
possibly required to improve our result. 


The dBZe-DV correlations are principal aspects rep- 
resentative of cloud microphysics for the Doppler radar 
measurement and simulation. The variables are influ- 
enced by certain relevant physical factors, as in Eqs. (1) 
and (8). Figures 11, 12, and 13 show the dBZe-DV re- 
lationships in the form of histograms and line graphs of the 
averaged DV with respect to the dBZe range of 1 dB. The 
observations and simulations share a common charac- 
teristic: the DV decreases (negative velocities are down- 
ward) with dBZe under a height of 4 km, whereas the DV 





2580 


JOURNAL OF THE ATMOSPHERIC SCIENCES 


Volume 69 



Radarsig [dBZe] 
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Fig. 13. dBZe-averaged DV (a) under heights of 4 km and (b) over 
heights of 5 km. Solid, dash-dotted, dash-double dotted, and dotted 
lines denote distributions by the radar measurement, the bin (con- 
trol and rimed snow) simulations, and the bulk model simulation, 
respectively. 

is nearly constant with respect to dBZe over 5 km. This 
feature is considered to be due to the differences of the 
uniformity of particle size and the dependence of the 
particle fall velocity on particle size between the ice-phase 
and liquid-phase clouds. A unimodal distribution in the 
histograms is observed over 5 km (Fig. 12). The bin model 
with snow fall velocity correction is the best simulation 
with respect to the occurrence and the averaged DV (Fig. 
13b). Under 4 km, a rough bimodal distribution centered 


at approximately (-5dB, -1 m s” 1 ) and (5 dB, -2.5 
ms -1 ) is noticed in the histogram plot of the observed 
data (Fig. 11). This distribution is organized by large 
raindrops included in the precipitation from the ice clouds 
and relatively small cloud droplets in the shallow water 
clouds. In contrast, the histograms of the simulations 
highlight only approximately one mode, as observed in the 
comparison of the DV vertical histograms (Fig. 8). In Fig. 
11, the profiles of the two bin model simulations show only 
the mode near (—5 dB, — 1 m s -1 ) in the observed profile; 
this result indicates a partial insufficiency of precipitation 
from the ice clouds also. The bulk simulation profile shows 
only another mode near (5 dB, -2.5 ms -1 ); the mode by 
shallow water clouds cannot be reproduced probably be- 
cause of slight dominance of the cloud water class in the 
bulk model simulation. 

b. CloudSat spaceborne CPR (C2873 and C2917) 

Comparing the simulations for the models and the 
CloudSat CPR observations may introduce another view- 
point to the issue. Ice clouds can be measured without at- 
tenuations by gas and hydrometeors under the target, 
which subsequently has a large influence on the calculation 
of the Ze of the ice clouds in the case of the shipborne 
observation. The dBZe profiles were illustrated for C2873 
and C2917 in the form of the latitude-height cross sections 
(LHCSs) (Figs. 14 and 15, respectively). The freezing level 
was approximately located at a height of 1-2 km in both 
cases. Figure 16 shows line charts for the vertical profile of 
the horizontally averaged dBZe with sampling from —40 to 
20 dB. The overestimations of dBZe in the ice cloud layers 
are highlighted, as in the case of M2223. 

For C2873, significant overestimations of PIA due to 
hydrometeors are highlighted, especially over the latitude 
from 40° to 44°N (Fig. 17). Ice clouds are a considerable 
cause of the overestimations because the simulated PIA 
from the top to a height of 2.5 km already exceeds the total 
PIA of the CloudSat product (not figured). The result 
demonstrates that the overestimations of dBZe and the 
attenuation are due to the characteristics of the simulated 
ice clouds themselves. Figures 14b and 14d show that part 
of Ze is less than —40 dB in the lower layer under 2 km. 
This feature is also attributed to strong attenuation in the 
ice clouds, although precipitation was actually simulated 
there. The adjustment in the additional run of the bin 
microphysics with the rimed snow is able to reduce the 
overestimations of PIA and Ze, leading to agreement with 
the observed data (Figs. 14c and 16a) because the vertical 
distribution of the dominant hydrometeor types (Fig. 18a) 
is similar to that in the M2223 case (Fig. 6a): Snow and ice 
crystals are dominant in the ice clouds of the bin model 
simulation, whereas only snow is dominant in the bulk 
model simulation. 
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d) Bulk model simulation 
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Fig. 14. Latitude-height cross sections of equivalent radar reflectivity factors in the case of 
nimbostratus on 11 Nov 2006 (a) measured by the CloudSat CPR, and calculated using the 
radar simulator applied to the outputs of (b) the bin (control), (c) the bin (rimed snow), and (d) 
the bulk model simulation. 


In the case of C2917, shallow convective clouds, in- 
cluding more supercooled water, were simulated. Con- 
sequently, hail and graupel were present in relatively 
high percentages compared with the other cases in both 
simulations (Figs. 18c, d). As a result, the adjustment of 
the bin microphysics with the riming did not significantly 
improve the discrepancy between the observed data and 
the simulation (Fig. 16b) because these modulations are 
mostly effective in changing the microphysical struc- 
tures of snow and ice crystals. In addition, an error in the 
simulations of the macrophysical structures of clouds 


may cause the overestimation in this case. The cloud-top 
heights were simulated to be over 4 km in the northern 
part of the latitude from 38.5°N, which is not in agree- 
ment with the observed data. The result suggests stronger 
convection, which causes increases in IWC, mean particle 
size, and, subsequently, dBZe. 

4. Conclusions 

In this paper, we conducted NWP simulations using 
JMA-NHM with bin and bulk microphysics for three 
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Fig. 15. As in Fig. 14, but for the case of convective clouds on 14 Nov 2006. 


cases to be compared with shipborne or spaceborne 
cloud radar observations. Direct comparisons were 
performed with a radar product simulator applied to the 
output of the NWP simulations. The results documented 
some important characteristics and problems of micro- 
physics simulation in the bin and bulk models. The bin 
model can generally provide better Ze and DV simula- 
tions than the bulk model. The modification of the ter- 
minal fall velocities of snowflakes equalized to those of 
hail in all size bins is considerably effective for improving 
the result. The result suggests that there is something 
wrong in the mass-size and size-terminal fall velocity 
relations of snowflakes in the present bin microphysics. 


The consequence is consistent with the result of the study 
using the same bin microphysics core (Fan et al. 2009). 
Insufficient modeling of the riming process is a possible 
source of the problem, although a significant amount of 
supercooled water is not predicted in the simulation. 
Advanced microphysical models explicitly calculating 
a rimed fraction of snow have the potential to provide 
a better simulation. A poor modeling of snow aggrega- 
tion is one of the possibilities but the required increase in 
fall velocity is considered not to be achieved using a dif- 
ferent shape assumption of snow aggregates without 
riming (Pruppacher and Klett 1997). Another possibility 
is an underestimation of increase in the terminal fall 
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Fig. 16. Vertical distributions of the horizontally averaged dBZe 
sampling the range from -40 to 20 dB in (a) Fig. 14 and (b) Fig. 15. 


velocity aloft by decrease in air pressure or density. In 
Lin et al. (1983), the square root factor in the equation 
is derived from the study of Foote and du Toit (1969) 
about the terminal velocity of raindrops aloft. It is un- 
certain whether the exponents of the factor for liquid- 
phase and ice-phase particles are the same, since their 
drag coefficients are different. In addition, there still 
remains a lack of understanding how temporal or spatial 
mismatches between measurement and simulation affect 
the result of comparison, although the effect can be 
relaxed in a quantitative analysis. An insufficiency of 
precipitation from ice clouds is observed in the one case, 
and the problem should be investigated in sampling 
many cases in future work. 

We provide some suggestions for improving the bulk 
microphysics that are based upon comparisons with the 
observed data and the bin model result. First, one of the 
problems in Ze simulation with the bulk microphysics is 
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Fig. 17. Horizontal distributions of the path integration attenua- 
tion (PIA) in the case of nimbostratus on 11 Nov 2006. 
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a lack of countered spread in the histogram for ice 
clouds. Bodas-Salcedo et al. (2008) suggested that it is 
necessary to predict the number concentration in a 
double-moment framework. However, at least in our 
cases, a lack of dBZe spread of ice clouds may be re- 
solved if different classes coexist adequately because 
snow is overdominant in the bulk model simulation 
compared with the bin model result. The overdominance 
of snow is considered to be due to the large depositional 
growth to snow (Eito and Aonashi 2009). Second, the 
vertical structure of the DV in the bulk model simulation 
is clearly different from those of the observed data and 
the bin model simulation. The fall velocity in each cloud 
class is dependent on the mixing ratio and air density. A 
decrease in DV with altitude due to changes in air density 
should be offset by increases due to changes in particle 
size; however, the increases do not function in the bulk 
model simulation. The effect should be introduced in the 
calculation of the fall velocity, even in single -moment 
bulk microphysics. A double-moment bulk microphysics 
can calculate the fall velocity by taking into account 
changes in the mean particle size, and hence the problem 
is also avoidable. 

An additional approach using multiwavelength remote 
sensing is useful for further study. Lidar observations 
are often accompanied by W-band radar observations, 
such as CALIPSO with CloudSat. A spaceborne lidar is 
ideal for detecting thin cirrus; a fraction of up to 70% of 
upper clouds detected by CloudSat can be measured by 
the lidar. Measuring very thick clouds using the lidar is 
difficult because they cannot be penetrated by the lasers. 
In such regions, the combination of X-band, C-band, or 
Ka-band radar with W-band could offer a potential 
solution. X-band, C-band, or K-band frequencies are 
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Fig. 18. Vertical distribution of the horizontally averaged mass ratios of each hydrometeor category to the total 
water contents in the CloudSat cases. (a),(c) Droplets (solid line), the total of three types of ice crystals (dash-dotted), 
snow (dotted), and the total of graupel and hail (dash-double dotted) in the bin simulation; (b),(d) the total of cloud 
water and rain (solid line), cloud ice (dash-dotted), snow (dotted), and graupel (dash-double dotted) in the bulk 
model simulation. 


sometimes used together with the W-band frequency; 
TRMM PR and CloudSat for global data constitute an 
example (Masunaga et al. 2008). K-band multiwave- 
length radars or radiometers are utilized to observe 
clouds on site (e.g., Illingworth et al. 2007) and are 
currently being prepared in the Global Precipitation 
Measurement (GPM) mission as a next-generation 
satellite-based observation. In addition, collaboration 
between radar remote sensing and on-site measure- 
ment is necessary to verify the microphysics. This kind 


of measurement campaign is often performed for the 
purpose of on-site validation of instruments for remote 
sensing. 
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